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Abstract 

This  research  had  as  its  primary  objective  the  development  of  rigorous  physical 
descriptions  and  corresponding  numerical  strategies  for  the  treatment  of  dynamically 
contacting  solid  and  structural  entities,  and  produced  useful  new  algorithms  for  de¬ 
scription  of  dynamic  frictional  response,  coupled  thermomechanical  interfacial  behav¬ 
iors,  elastic  and  inelastic  impact  interaction,  and  tribological  complexity.  The  unifying 
goal  of  this  work  was  to  assure  that  problems  involving  transient  impacts  between 
solids  could  be  numerically  solved  with  the  same  level  of  robustness,  accuracy,  and 
physical  insight  as  is  currently  brought  to  bear  in  other  areas  of  nonlinear  mechanics. 
The  results  obtained  in  this  work  have  importance  well  beyond  contact  mechanics, 
extending  also  to  other  applications  where  computational  models  of  neighboring,  but 
distinct,  domains  must  be  coupled  together  in  a  numerically  accurate  manner. 
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1  Objectives  and  Accomplishments  of  this  Research 

This  research  had  as  its  overarching  goal  the  development  of  contact/impact  algorithms 
for  computational  nonlinear  solid  mechanics,  such  that  finite  element  simulation  of  interfa¬ 
cial  phenomena  could  be  accomplished  with  the  same  degree  of  robustness,  accuracy,  and 
reliability  currently  present  in  algorithmic  treatments  of  such  features  as  inelasticity,  large 
deformation  response,  and  dynamics.  The  effort  has  relevance  not  only  to  contact  mechanics, 
but  also  to  other  applications  (such  as  phase  transformation,  surface  adhesion  and  debond¬ 
ing,  fragmentation  and  fracture,  and  structural  composite  response)  where  distinct  domains 
must  be  coupled  together  in  a  numerically  accurate  manner. 

Specific  achievements  of  this  research,  fitting  under  the  aforementioned  overall  objective, 
included: 

•  Reliable  treatments  of  dynamic  frictional  response,  including  new  interface  constitutive 
laws  that  describe  slip  instability,  as  well  as  corresponding  numerical  algorithms  that 
treat  such  models  in  a  stable  manner 

•  Improved  numerical  description  of  coupled  thermomechanical  interfacial  behavior,  treat¬ 
ing  in  particular  frictional  heating  and  thermal  softening  associated  with  friction  and 
developing  both  monolithic  and  partitioned  solution  strategies  for  such  problems; 

•  New  finite  element  algorithms  for  impact  interaction,  considering  the  energy-momentum 
paradigm  as  a  new  framework  in  which  unconditional  stability  of  impact  calculations 
could  be  assured  and  in  which  partitioning  of  energy  and  momentum  during  impact 
events  (both  elastic  and  inelastic)  was  accurately  predicted;  and 

•  New  methods  for  spatial  discretization  of  contact  problems,  developing  new  means  by 
which  stresses  on  interfaces  could  be  accurately  approximated  and  by  which  dissimilar 
domains  can  be  reliably  joined  numerically. 

Subsequent  sections  wall  describe  advances  made  in  each  of  these  areas  in  somewhat  more 
detail. 

2  Technological  Issues  Considered  in  this  Project 

2.1  Computational  Treatment  of  Stick-Slip  Instabilities  in  Non¬ 
linear  Mechanics 

It  has  been  known  for  centuries  that  under  conditions  of  constant  pressure  across  an  inter¬ 
face,  the  frictional  force  required  to  sustain  kinetic  frictional  sliding  is  generally  lower  than 
that  required  to  initiate  it  (see,  e.g.,  [Eul48]).  Indeed,  it  is  also  often  observed  that  the 
apparent  coefficient  of  friction  prevailing  on  an  interface  is  dependent  on  sliding  velocity,  al¬ 
though  for  many  metals,  over  at  least  certain  ranges  of  sliding  velocity,  this  variation  occurs 
slowly  enough  so  that  the  “kinetic”  coefficient  of  friction  might  be  taken  as  a  constant,  in 
general  lower  than  the  “static”  coefficient  ([Rab65]).  In  general,  the  capability  to  introduce 
functional  rate  dependence  of  frictional  stress  is  also  desirable  in  many  applications. 
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Figure  1:  Tube  drawing  simulation:  (a)  schematic  representation  of  the  drawing  procedure 
and  initial  configuration  of  the  mesh;  (b)  final  configuration  of  the  mesh. 


Complicating  the  picture  somewhat,  however,  is  the  experimentally  observed  fact  that  a 
simple  one-to-one  relation  between  friction  coefficient  and  instantaneous  sliding  rate  is  not 
appropriate  for  many  interface  combinations.  In  particular,  as  examples  presented  in  [OM85] 
show,  the  apparent  friction  coefficient  in  a  dynamic  setting  is  often  not  a  reversible  function 
of  the  sliding  velocity,  but  displays  hysteretic  behavior  instead. 

Our  work  has  considered  the  use  of  rate  and  state  dependent  friction  laws,  of  the  sort 
pioneered  by  [Die79]  and  [Rui83]  in  the  geophysics  community,  for  engineering  applications 
of  nonlinear  mechanics.  The  models  developed  here  have  the  following  two  characteristics: 

1.  the  apparent  coefficient  of  friction  at  a  given  steady  state  velocity  is  a  single- valued 
(generally  decreasing)  function  of  that  velocity  through  the  state  variable  definition; 
and 

2.  responses  to  more  sudden  changes  in  slip  velocity  occur  on  a  length  scale  related  to 
typical  asperity  spacings  on  the  surface,  with  the  relationship  between  friction  stress 
and  slip  rate  being  of  a  viscous  nature  over  short  time  scales 

A  typical  application  of  our  methodology  is  shown  in  Figures  1  and  2,  which  collectively 
depict  a  tube  drawing  example.  Of  particular  note  is  the  inability  of  a  traditional  Coulomb 
description  to  describe  intermittent  “stick-slip”  frictional  behavior,  which  manifests  itself  in 
this  case  in  experimentally  observed  chatter  of  the  drawn  part  against  the  machine  tools. 
The  rate  and  state  dependent  friction  model  developed  in  this  research  reproduces  the  ex¬ 
perimentally  observed  behavior,  and  in  the  presence  of  large  deformations  and  significant 
inelasticity. 


4 


50 


Time(sec) 


(a) 


45 


25  1 - * - * - - - * - 1 

2.0  2.2  2.4  2.6 

Time(sec) 

(b) 


Figure  2:  Tube  drawing  simulation:  (a)  computed  force  on  the  die  during  early  deformation, 
with  experimental  data  of  [Coo81]  shown  in  inset;  (b)  computed  drawing  force  on  the  tube, 
after  some  length  has  already  been  drawn. 

2.2  Robust  Algorithmic  Strategies  for  Coupled  Thermomechani¬ 
cal  Friction  on  Interfaces 

Another  product  of  this  research  was  been  the  extension  of  computational  methods  to  encom¬ 
pass  thermomechanical  effects  on  interfaces,  such  that  frictional  heating,  thermal  frictional 
softening,  and  coupling  effects  are  explicitly  considered  by  the  formulation.  A  full  treat¬ 
ment  of  both  thermal  and  mechanical  effects  was  included  in  our  modeling  approach,  such 
that  balance  of  momentum  and  balance  of  energy  were  enforced  throughout  the  contacting 
continua. 

A  key  factor  in  the  approach  followed  in  this  research  was  the  development  of  discrete 
numerical  representations  obeying  underlying  thermodynamical  restrictions.  In  particular, 
a  priori  stability  estimates  were  used  to  guide  the  development  of  integrators  for  the  fric¬ 
tional  equations  of  evolution,  and  to  extend  such  estimates  to  the  thermally  coupled  case 
so  that  unconditionally  stable  staggered  algorithms  for  thermomechanical  friction  problems 
resulted.  The  overall  algorithmic  approach  was  similar  to  that  utilized  by  [AS93]  in  their 
study  of  thermoplasticity  in  bulk  materials,  but  required  development  of  the  appropriate 
interface  kinematics  and  thermodynamical  restrictions  to  enable  new  stability  estimates  to 
be  constructed  and  numerically  approximated  (see  [Lau99]).  An  example  benefiting  from  the 
technology  developed  is  shown  in  Figures  3  and  4,  which  depict  a  finite  element  simulation 
of  a  Howitzer  shell  firing.  Here,  the  frictional  heating  and  subsequent  thermal  softening  of 
the  shell/barrel  interface  is  a  key  factor  in  the  performance  of  the  system  during  firing. 
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Figure  3:  Frictional  heating  in  a  155  mm  Howitzer:  (a)  schematic  representation  of  projectile 
and  dimensions  of  the  band;  (b)  undeformed  initial  mesh  for  the  band  and  tube  (upper)  and 
deformed  configuration  and  temperature  map  after  50  mm  of  sliding. 
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Figure  4:  Projectile-barrel  problem:  comparison  with  available  experimental  data,  (a)  Evo¬ 
lution  of  the  average  apparent  coefficient  of  friction  as  the  projectile  advances  along  the 
tube:  finite  element  simulation  and  experimental  data  compiled  from  [Mon76];  (b)  average 
temperature  evolution  at  the  surface  of  the  band  as  sliding  proceeds  along  the  tube,  showing 
also  the  temperature  distribution  for  the  tube  after  500  mm  of  sliding;  (c)  average  apparent 
coefficient  of  friction  as  a  function  of  temperature. 
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2.3  New  Energy  Consistent  Algorithms  for  Inelastic  Impact  Sim¬ 
ulation 

The  new  paradigm  we  have  developed  for  inelastic  impact  simulation  has  its  origins  in  the 
energy-momentum  method  for  elastodynamics,  as  originally  proposed  by  Simo  and  Tarnow 
[ST92].  As  originally  conceived,  the  method  was  intended  to  ensure  unconditional  stability  in 
the  fully  nonlinear  regime,  while  also  providing  for  exact  conservation  of  linear  and  angular 
momentum  for  autonomous  dynamical  systems.  In  the  work  done  under  this  contract,  the 
paradigm  has  been  extended  not  only  to  provide  for  unconditionally  stable  treatment  of  im¬ 
pact  interactions,  but  also  to  allow  for  inelastic  dissipation  (giving  rise  to  energy  consistency, 
rather  than  conservation)  both  on  contact  interfaces  and  in  the  bulk  continua  (see  [ML02], 
Furthermore,  introduction  of  temporal  discontinuities  in  the  velocity-displacement  updates 
utilized  on  contact  interfaces,  in  a  manner  consistent  with  the  physics  being  simulated  and 
the  conservation  laws  being  enforced,  has  produced  a  new  velocity  update  algorithm  [LL01] 
which  provides  arbitrarily  accurate  enforcement  of  contact  constraints,  giving  an  unprece¬ 
dented  combination  of  stability  and  accuracy. 

In  impact  mechanics,  each  local  contact  event  introduces  a  jump  in  the  contact  velocities, 
with  the  magnitude  of  the  jump  being  related  to  the  material  properties  of  the  contacting 
entities  (in  an  elastic  analysis,  these  would  be  the  impedances  of  the  bodies  in  question). 
Our  work  has  shown  that  algorithmic  counterparts  of  these  jumps  can  be  introduced  into 
the  velocity-displacement  relationship  on  interfaces,  with  the  magnitude  of  the  jumps  being 
found  from  the  global  conservation  conditions.  The  algorithm  is  applicable  to  any  given 
constitutive  description,  and  is  designed  for  fully  nonlinear  applications  in  large  deformations. 
It  can  enforce  the  impenetrability  conditions  associated  with  contact  arbitrarily  closely,  while 
also  retaining  all  conservation  properties.  An  example  is  given  in  the  two  ring  impact  example 
in  Figure  5,  which  depicts  also  the  accurate  enforcement  of  the  impenetrability  conditions  in 
a  typical  time  step  (inset).  The  evolution  of  the  total  system  energy  and  angular  momentum 
is  shown  in  Figure  6. 

Another  key  achievement  pertaining  to  this  portion  of  the  project  was  the  development  of 
an  energy  consistent  algorithm  for  finite  strain  plasticity,  such  that  bulk  inelasticity  could  be 
simulated  within  the  energy  momentum  framework.  Indeed,  our  previous  work  on  frictional 
contact  algorithms  had  shown  that  dissipative  phenomena  could  be  treated  within  the  energy 
momentum  framework  by  forcing  a  discrete  version  of  an  a  priori  stability  estimate  to  be 
algorithmically  obeyed.  In  this  research,  we  succeeded  in  implementing  this  idea  for  two 
classes  of  inelastic  solids:  an  additively  decomposed  hyperelastic-plastic  constitutive  law 
based  on  a  plastic  Green’s  strain,  and  a  model  based  on  the  multiplicative  decomposition  of 
the  deformation  gradient  into  elastic  and  plastic  components.  An  inelastic  impact  simulation 
benefiting  from  this  methodology  is  shown  in  Figures  7  and  8;  here,  it  is  seen  that  a  linearly 
stable  time  integrator  can  produce  unstable  results  in  the  fully  nonlinear  regime  even  in  the 
presence  of  potentially  significant  inelastic  dissipation.  The  energy  consistent  algorithms 
developed  in  this  work  are  seen  to  eliminate  this  possibility,  producing  dissipation  in  a 
manner  consistent  with  the  material  and  interface  laws  in  use. 
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Figure  5:  Two  ring  impact  problem,  dynamic  configurations,  t  —  [0, 2,4,  ...,20].  The  inset 
shows  a  zoom  of  the  contact  interface  at  a  typical  time  step.  The  velocity  update  algorithm 
is  seen  to  produce  stable  results  and  to  facilitate  accurate  constraint  enforcement. 
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Figure  6:  Two  ring  impact  problem:  (a)  energy  evolution,  (b)  angular  momentum  evolution. 
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Figure  7:  Deformed  sequence  corresponding  to  an  elastoplastic  ring  impact  calculation,  using 
the  energy  consistent  algorithm  developed  in  this  research. 
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Figure  8:  Energy  evolution  in  the  elastoplastic  ring  impact  problem,  using  energy  consistent 
and  Newmark  algorithms:  (a)  the  whole  history  of  evolution;  (b)  enlarged  plot  showing 
energy  evolution  during  impact. 
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Figure  9:  Discretization,  cylinder  on  cylinder  contact  problem. 


2.4  New  Spatial  Discretization  Strategies  for  Contact  Interaction 

The  mortar  element  method,  previously  developed  for  domain  decomposition  problems,  was 
extended  in  this  project  to  provide  a  new  technique  for  the  spatially  discretization  of  contact 
interactions.  The  method  relies  on  an  intermediate  contact  surface,  over  which  all  contact 
variables  (kinematic  and  traction)  are  defined.  These  contact  variables  are  related  via  a 
variational  projection  to  their  counterparts  on  the  two  contacting  bodies,  such  that  the 
contact  constraints  (enforced  over  the  intermediate  surface)  are  related  in  a  nonlocal  manner 
to  the  displacement  and  force  quantities  on  the  surfaces  of  the  contacting  bodies.  This 
procedure  results  in  an  accurate  and  variationally  consistent  algorithm  to  recover  contact 
pressures  on  contact  surfaces,  and  ensures  passage  of  contact  patch  tests  in  general.  A 
frictional  example  depicting  the  effectiveness  of  the  procedure  is  shown  in  Figures  9  and  10, 
where  a  frictional  Hertzian  example  is  considered.  As  can  been  seen,  the  contact  tractions 
are  very  accurately  represented  using  this  new  procedure,  even  though  the  meshes  of  the 
contacting  bodies  are  nonconforming  throughout  the  contact  patch. 

Another  technological  issue  associated  with  spatial  discretization  of  contact  interaction 
is  seen  when  tangential  relative  motion  occurs  between  two  contacting  boundaries.  In  this 
instance,  sliding  must  take  place  in  general  across  nonsmooth  regions  created  by  the  C°  finite 
element  discretization,  which  gives  rise  to  a  lack  of  numerical  robustness  and/or  oscillations 
in  the  contact  tractions.  For  any  reasonable  description  of  tribological  complexity  (e.g., 
wear,  lubrication)  to  take  place,  these  issues  of  spatial  discretization  should  be  resolved. 
This  research  succeeded  in  developing  a  new  Hermitian  smoothing  method  for  contact  inter¬ 
actions,  relying  on  an  underlying  variational  principle  to  consider  alternative  (geometrically 
smoothed)  representations  of  contact  geometry.  The  result  is  greatly  enhanced  numerical 
robustness,  as  well  as  a  significant  decrease  in  observed  numerical  oscillations.  An  example 
of  the  application  of  such  ideas  is  given  in  Figures  11  and  12,  which  depict  a  drawbead 
clamping  operation  such  as  is  commonly  encountered  in  metal  forming.  As  can  be  seen,  the 


Figure  10:  Computed  mortar  surface  nodal  tractions  for  frictional  cylinder  on  cylinder  con¬ 
tact  problem. 


contact  smoothing  procedure  gives  rise  to  a  significant  improvement  in  the  numerical  results 
obtained,  giving  an  improvement  also  by  about  a  factor  of  four  in  computational  time. 


3  Statement  of  Relevance  to  the  Navy 

Among  the  Navy  applications  requiring  good  interfacial  descriptions  are:  weapons  effect 
calculations,  durability  studies,  predictions  of  various  machining  and  materials  processing 
operations,  characterization  of  nonlinear  vibrational  behavior  of  various  onboard  systems, 
and  even  description  of  solid/fluid  interaction,  where  the  solid/fluid  interface  may  be  thought 
of  as  a  contact  interface.  Without  numerical  algorithms  in  which  temporal  and  spatial  accu¬ 
racy  are  assured,  the  predictive  power  of  such  calculations  remains  an  open  question.  This 
research  has  made  significant  algorithmic  and  theoretical  contributions  which  significantly 
enhance  our  ability  to  perform  these  types  of  calculations.  In  this  sense,  it  seems  that  the 
goals  of  the  project  as  stated  in  the  original  proposal  have  been  achieved. 
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Figure  11:  Undeformed  configuration  of  the  sheet-drawbead  system. 


Figure  12:  Force-displacement  plots  for  the  smoothed  and  unsmoothed  drawbead  problems 
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